|
|
(*~\Глагол\Отделы\Числа~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
(**) ОТДЕЛ Функц;
(*============================================================================*
* НАЗНАЧЕНИЕ: специальные функции
*----------------------------------------------------------------------------*
* ПРИМЕЧАНИЯ: в особых случаях в <Матем.код> записывается # 0 значение.
* Источники, ссылки, библиография
*============================================================================*)
ИСПОЛЬЗУЕТ
Матем;
ВИД
Вещ = Матем.Вещ;
ТаблКоэф = РЯД 10 ИЗ Вещ;
ПОСТ
ТАБЛФАКТ = 34; (* размер таблицы, хранящей факториалы небольших чисел *)
МАКСФАКТ = 170; (* наибольший аргумент для вычисления факториала *)
(* для НепГ и НепБета *)
BIG = 9.223372036854775808D18;
BIGINV = 1.084202172485504434007D-19;
(* для Г функций *)
МАКСГ = 171.624376956302D0;
МАКСЛНГ = 2.556348D+305;
ПЕР
STIR,STIRP,GAMS,GAMSN,GAMP,GAMQ,LNGAMP,LNGAMQ:ТаблКоэф;
INP0,INP1,INP2,INP3:ТаблКоэф;
INQ0,INQ1,INQ2,INQ3:ТаблКоэф;
ТаблФакт:РЯД ТАБЛФАКТ ИЗ Вещ; (* таблица факториалов небольших чисел *)
i:ЦЕЛ;
(*-оооооооооооооо- опережающие определения/заготовки задач -ооооооооооооооооо-*)
ЗАДАЧА^ Факториал-(x:Вещ):Вещ; (* факториал *)
ЗАДАЧА^ БиномКоэф-(x,y:Вещ):Вещ; (* биномиальные коэффициенты C(n,k) *)
ЗАДАЧА^ Г-(x:Вещ):Вещ; (* Гамма функция *)
ЗАДАЧА^ ЗнакГ-(x:Вещ):ЦЕЛ; (* Знак Гамма функции - Г(x) *)
ЗАДАЧА^ lnГ-(x:Вещ):Вещ; (* ln(|Г(x)|) *)
ЗАДАЧА^ НепГ-(a,x:Вещ):Вещ; (* неполная Гамма функция *)
ЗАДАЧА^ ДопГ-(a,x:Вещ):Вещ; (* дополнение НепГ *)
ЗАДАЧА^ Бета-(x,y:Вещ):Вещ; (* Бета функция *)
ЗАДАЧА^ НепБета-(a,b,x:Вещ):Вещ; (* неполная Бета функция *)
ЗАДАЧА^ erf-(x:Вещ):Вещ; (* функция ошибок *)
ЗАДАЧА^ erfc-(x:Вещ):Вещ; (* дополнительная функция ошибок *)
(*----------------------------------------------------------------------------*
* Биномиальное распределение с вероятностью < p > и числом испытаний < n >
*----------------------------------------------------------------------------*)
ЗАДАЧА^ БиномВ-(n:ЦЕЛ; p:Вещ; k:ЦЕЛ):Вещ; (* p(x = k) *)
ЗАДАЧА^ БиномП-(n:ЦЕЛ; p:Вещ; k:ЦЕЛ):Вещ; (* p(x <= k) *)
(*----------------------------------------------------------------------------*
* Распределение Пуассона со средним значением < Mu >
*----------------------------------------------------------------------------*)
ЗАДАЧА^ ПуассонВ-(Mu:Вещ; k:ЦЕЛ):Вещ; (* p(x = k) *)
ЗАДАЧА^ ПуассонП-(Mu:Вещ; k:ЦЕЛ):Вещ; (* p(x <= k) *)
(*----------------------------------------------------------------------------*
* Стандартный нормальный закон распределения
*----------------------------------------------------------------------------*)
ЗАДАЧА^ НормПл-(x:Вещ):Вещ; (* Плотность распределения *)
ЗАДАЧА^ НормП-(x:Вещ):Вещ; (* p( U <= x) *)
ЗАДАЧА^ НормВ-(x:Вещ):Вещ; (* p(|U|>=|x|) *)
ЗАДАЧА^ ОбрНорм-(p:Вещ):Вещ; (* Возвращает такое , что p(U <= x) = < p > *)
(*----------------------------------------------------------------------------*
* Распределение Стьюдента с < Nu > степенями свободы
*----------------------------------------------------------------------------*)
ЗАДАЧА^ СтьюдентПл-(Nu:ЦЕЛ; x:Вещ):Вещ; (* Плотность распределения *)
ЗАДАЧА^ СтьюдентП-(Nu:ЦЕЛ; x:Вещ):Вещ; (* p( t <= x) *)
ЗАДАЧА^ СтьюдентВ-(Nu:ЦЕЛ; x:Вещ):Вещ; (* p(|t|>=|x|) *)
(*----------------------------------------------------------------------------*
* Распределение Хи квадрат с < Nu > степенями свободы
*----------------------------------------------------------------------------*)
ЗАДАЧА^ Хи2Пл-(Nu:ЦЕЛ; x:Вещ):Вещ; (* Плотность распределения Хи2 *)
ЗАДАЧА^ Хи2П-(Nu:ЦЕЛ; x:Вещ):Вещ; (* p(Khi2 <= x) *)
ЗАДАЧА^ Хи2В-(Nu:ЦЕЛ; x:Вещ):Вещ; (* p(Khi2 >= x) *)
(*----------------------------------------------------------------------------*
* Распределение Фишера с < Nu1 > и < Nu2 > степенями свободы
*----------------------------------------------------------------------------*)
ЗАДАЧА^ ФишерПл-(Nu1,Nu2:ЦЕЛ; x:Вещ):Вещ; (* Плотность распределения *)
ЗАДАЧА^ ФишерП-(Nu1,Nu2:ЦЕЛ; x:Вещ):Вещ; (* p(f <= x) *)
ЗАДАЧА^ ФишерВ-(Nu1,Nu2:ЦЕЛ; x:Вещ):Вещ; (* p(f >= x) *)
(*----------------------------------------------------------------------------*
* Экспоненциальное распределение
*----------------------------------------------------------------------------*)
ЗАДАЧА^ ЭкспПл-(a,x:Вещ):Вещ; (* Плотность распределения *)
ЗАДАЧА^ ЭкспП-(a,x:Вещ):Вещ; (* p( <= x) *)
(*----------------------------------------------------------------------------*
* Бета распределение
*----------------------------------------------------------------------------*)
ЗАДАЧА^ БетаПл-(a,b,x:Вещ):Вещ; (* Плотность распределения *)
ЗАДАЧА^ БетаП-(a,b,x:Вещ):Вещ; (* p( <= x) *)
(*----------------------------------------------------------------------------*
* Г-распределение
*----------------------------------------------------------------------------*)
ЗАДАЧА^ ГПл-(a,b,x:Вещ):Вещ; (* Плотность распределения *)
ЗАДАЧА^ ГП-(a,b,x:Вещ):Вещ; (* p( <= x) *)
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Горнер0(x:Вещ; a-:ТаблКоэф; n:ЦЕЛ):Вещ;
(*============================================================================*
(* Цель: вычислить многочлен степени < n >
*
* y = a0xn + a1xn-1 +...+ an-1x1 + an
*
* по схеме Горнера.
*============================================================================*)
ПЕР
i:ЦЕЛ;
ответ:Вещ;
УКАЗ
ответ:=a[0];
ОТ i:=1 ДО n ВЫП
ответ:=ответ*x+a[i]
КОН;
ВОЗВРАТ ответ
КОН Горнер0;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Горнер1(x:Вещ; a-:ТаблКоэф; n:ЦЕЛ):Вещ;
(*============================================================================*
(* Цель: вычислить многочлен степени
*
* y = xn + a0xn-1 +...+ an-2x1 + an-1
*
* по схеме Горнера.
*============================================================================*)
ПЕР
ответ:Вещ;
i:ЦЕЛ;
УКАЗ
ответ:=x+a[0];
ОТ i:=1 ДО n-1 ВЫП
ответ:=ответ*x+a[i]
КОН;
ВОЗВРАТ ответ
КОН Горнер1;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ЗнакГ-(x:Вещ):ЦЕЛ;
УКАЗ
ЕСЛИ (x > 0) ИЛИ НЕ ЧЕТ(ВШИРЦЕЛ(МОДУЛЬ(x))) ТО
ВОЗВРАТ 1
КОН;
ВОЗВРАТ -1
КОН ЗнакГ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Stirf(x:Вещ):Вещ;
(* Формула Стирлинга для Гамма функции:
* Г(x) = Sqrt(2*Pi) * x^((x-0.5)*exp(-x)) * (1+(1/x)*p(1/x))
* где p(x) является многочленом *)
ПЕР
W,p:Вещ;
УКАЗ
W:=1/x;
ЕСЛИ x > 1024 ТО
p:=6.97281375836585777429D-5*W + 7.84039221720066627474D-4;
p:=p*W - 2.29472093621399176955D-4;
p:=p*W - 2.68132716049382716049D-3;
p:=p*W + 3.47222222222222222222D-3;
p:=p*W + 8.33333333333333333333D-2;
ИНАЧЕ
p:=Горнер0(W,STIR,8)
КОН;
ВОЗВРАТ Матем.КВКОР2ПИ*Матем.exp((x-0.5)*Матем.ln(x)-x)*(1+W*p)
КОН Stirf;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА GamSmall(x1,z:Вещ):Вещ;
(* Гамма функция для малых аргументов *)
ПЕР
p:Вещ;
УКАЗ
ЕСЛИ x1 = 0 ТО
Матем.код:=Матем.ОСОБЕННОСТЬ;
ВОЗВРАТ Матем.МАКСВЕЩ
КОН;
ЕСЛИ x1 < 0 ТО
x1:=-x1;
p:=Горнер0(x1,GAMSN,8)
ИНАЧЕ
p:=Горнер0(x1,GAMS,8)
КОН;
ВОЗВРАТ z/(x1*p)
КОН GamSmall;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА StirfL(x:Вещ):Вещ;
(* Приближение ln(Г) формулой Стирлинга при x >= 13 *)
ПЕР
q,W:Вещ;
УКАЗ
q:=Матем.ln(x)*(x-0.5)-x;
q:=q+Матем.ЛНКВКОР2ПИ;
ЕСЛИ x > 1.E+10 ТО
ВОЗВРАТ q
ИНАЧЕ
W:=1/Матем.кв(x);
ВОЗВРАТ q+Горнер0(W,STIRP,6)/x
КОН
КОН StirfL;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Г-(x:Вещ):Вещ;
ПЕР
SgnGam,n:ЦЕЛ;
a,x1,z:Вещ;
УКАЗ
Матем.код:=0;
SgnGam:=ЗнакГ(x);
ЕСЛИ (x = 0) ИЛИ ((x < 0) И (ЦЕЛЧАСТЬ(x) = x)) ТО
Матем.код:=Матем.ОСОБЕННОСТЬ;
ВОЗВРАТ SgnGam*Матем.МАКСВЕЩ
КОН;
ЕСЛИ x > МАКСГ ТО
Матем.код:=Матем.ПЕРЕПОЛНЕНИЕ;
ВОЗВРАТ Матем.МАКСВЕЩ
КОН;
a:=МОДУЛЬ(x);
ЕСЛИ a > 13 ТО
ЕСЛИ x < 0 ТО
n:=УЗК(ВШИРЦЕЛ(a));
z:=a-n;
ЕСЛИ z > 0.5 ТО
n:=n+1;
z:=a-n
КОН;
z:=МОДУЛЬ(a*Матем.sin(Матем.ПИ*z))*Stirf(a);
ЕСЛИ z <= Матем.ПИ/Матем.МАКСВЕЩ ТО
Матем.код:=Матем.ПЕРЕПОЛНЕНИЕ;
ВОЗВРАТ SgnGam*Матем.МАКСВЕЩ
КОН;
z:=Матем.ПИ/z
ИНАЧЕ
z:=Stirf(x)
КОН;
ВОЗВРАТ SgnGam*z
ИНАЧЕ
z:=1;
x1:=x;
ПОКА x1 >= 3 ВЫП
x1:=x1-1;
z:=z*x1
КОН;
ПОКА x1 < -0.03125 ВЫП
z:=z/x1;
x1:=x1+1
КОН;
ЕСЛИ x1 <= 0.03125 ТО
ВОЗВРАТ GamSmall(x1,z)
ИНАЧЕ
ПОКА x1 < 2 ВЫП
z:=z/x1;
x1:=x1+1
КОН;
ЕСЛИ (x1 = 2) ИЛИ (x1 = 3) ТО
ВОЗВРАТ z
ИНАЧЕ
x1:=x1-2;
ВОЗВРАТ z*Горнер0(x1,GAMP,7)/Горнер0(x1,GAMQ,8)
КОН
КОН
КОН
КОН Г;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА lnГ-(x:Вещ):Вещ;
ПЕР
n:ЦЕЛ;
a,x1,z:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (x = 0) ИЛИ ((x < 0) И (ЦЕЛЧАСТЬ(x) = x)) ТО
Матем.код:=Матем.ОСОБЕННОСТЬ;
ВОЗВРАТ Матем.МАКСВЕЩ
КОН;
ЕСЛИ x > МАКСЛНГ ТО
Матем.код:=Матем.ПЕРЕПОЛНЕНИЕ;
ВОЗВРАТ Матем.МАКСВЕЩ
КОН;
a:=МОДУЛЬ(x);
ЕСЛИ a > 34 ТО
ЕСЛИ x < 0 ТО
n:=УЗК(ВШИРЦЕЛ(a));
z:=a-n;
ЕСЛИ z > 0.5 ТО
n:=n+1;
z:=n-a;
КОН;
z:=a*Матем.sin(Матем.ПИ*z);
ЕСЛИ z = 0 ТО
Матем.код:=Матем.ПЕРЕПОЛНЕНИЕ;
ВОЗВРАТ Матем.МАКСВЕЩ
КОН;
z:=Матем.ЛНПИ-Матем.ln(z)-StirfL(a);
ИНАЧЕ
z:=StirfL(x)
КОН;
ВОЗВРАТ z
АЕСЛИ x < 13 ТО
z:=1;
x1:=x;
ПОКА x1 >= 3 ВЫП
x1:=x1-1;
z:=z*x1
КОН;
ПОКА x1 < 2 ВЫП
ЕСЛИ МОДУЛЬ(x1) <= 0.03125 ТО
ВОЗВРАТ Матем.ln(МОДУЛЬ(GamSmall(x1,z)))
КОН;
z:=z/x1;
x1:=x1+1
КОН;
ЕСЛИ z < 0 ТО
z:=-z
КОН;
ЕСЛИ x1 = 2 ТО
ВОЗВРАТ Матем.ln(z)
ИНАЧЕ
x1:=x1-2;
ВОЗВРАТ x1*Горнер0(x1,LNGAMP,6)/Горнер1(x1,LNGAMQ,7)+Матем.ln(z)
КОН;
ИНАЧЕ
ВОЗВРАТ StirfL(x)
КОН
КОН lnГ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА НепГ-(a,x:Вещ):Вещ;
ПЕР
ответ,Ax,C,r:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (x <= 0) ИЛИ (a <= 0) ТО
ВОЗВРАТ 0
КОН;
ЕСЛИ (x > 1) И (x > a) ТО
ВОЗВРАТ 1-ДопГ(a,x)
КОН;
Ax:=a*Матем.ln(x)-x-lnГ(a);
ЕСЛИ Ax < Матем.МИНЛН ТО
Матем.код:=Матем.ОБНУЛЕНИЕ;
ВОЗВРАТ 0
КОН;
Ax:=Матем.exp(Ax);
(* степенной ряд *)
r:=a;
C:=1;
ответ:=1;
ПОВТОРЯТЬ
r:=r+1;
C:=C*x/r;
ответ:=ответ+C;
ДО C/ответ <= Матем.МАШЕПС;
ВОЗВРАТ ответ*Ax/a;
КОН НепГ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ДопГ-(a,x:Вещ):Вещ;
ПЕР
ответ,C,Yc,Ax,y,z,r,T,
Pk,Pkm1,Pkm2,Qk,Qkm1,Qkm2:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (x <= 0) ИЛИ (a <= 0) ТО
ВОЗВРАТ 1
КОН;
ЕСЛИ (x < 1) ИЛИ (x < a) ТО
ВОЗВРАТ 1-НепГ(a,x)
КОН;
Ax:=a*Матем.ln(x)-x-lnГ(a);
ЕСЛИ Ax < Матем.МИНЛН ТО
Матем.код:=Матем.ОБНУЛЕНИЕ;
ВОЗВРАТ 0
КОН;
Ax:=Матем.exp(Ax);
(* повторяющаяся дробь *)
y:=1-a;
z:=x+y+1;
C:=0;
Pkm2:=1;
Qkm2:=x;
Pkm1:=x+1;
Qkm1:=z*x;
ответ:=Pkm1/Qkm1;
ПОВТОРЯТЬ
C:=C+1;
y:=y+1;
z:=z+2;
Yc:=y*C;
Pk:=Pkm1*z-Pkm2*Yc;
Qk:=Qkm1*z-Qkm2*Yc;
ЕСЛИ Qk # 0 ТО
r:=Pk/Qk;
T:=МОДУЛЬ((ответ-r)/r);
ответ:=r
ИНАЧЕ
T:=1
КОН;
Pkm2:=Pkm1;
Pkm1:=Pk;
Qkm2:=Qkm1;
Qkm1:=Qk;
ЕСЛИ МОДУЛЬ(Pk) > BIG ТО
Pkm2:=Pkm2/BIG;
Pkm1:=Pkm1/BIG;
Qkm2:=Qkm2/BIG;
Qkm1:=Qkm1/BIG
КОН
ДО T <= Матем.МАШЕПС;
ВОЗВРАТ ответ*Ax
КОН ДопГ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Факториал-(x:Вещ):Вещ;
ПЕР
n:ЦЕЛ;
УКАЗ
n:=УЗК(ВШИРЦЕЛ(x));
Матем.код:=0;
ЕСЛИ n < 0 ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 1
КОН;
ЕСЛИ n > МАКСФАКТ ТО
Матем.код:=Матем.ПЕРЕПОЛНЕНИЕ;
ВОЗВРАТ Матем.МАКСВЕЩ
КОН;
ЕСЛИ n < ТАБЛФАКТ ТО
ВОЗВРАТ ТаблФакт[n]
ИНАЧЕ
ВОЗВРАТ Г(n+1)
КОН
КОН Факториал;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА БиномКоэф-(x,y:Вещ):Вещ;
ПЕР
n,k,i,n1:ЦЕЛ;
произв:Вещ;
УКАЗ
n:=УЗК(ВШИРЦЕЛ(x));
k:=УЗК(ВШИРЦЕЛ(y));
Матем.код:=0;
ЕСЛИ k < 0 ТО
ВОЗВРАТ 0
КОН;
ЕСЛИ (k = 0) ИЛИ (k = n) ТО
ВОЗВРАТ 1
КОН;
ЕСЛИ (k = 1) ИЛИ (k = n-1) ТО
ВОЗВРАТ n
КОН;
ЕСЛИ k > n-k ТО k:=n-k КОН;
n1:=n+1;
произв:=n;
ОТ i:=2 ДО k ВЫП
произв:=произв*((n1-i)/i)
КОН;
ВОЗВРАТ ВШИРЦЕЛ(произв)
КОН БиномКоэф;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Бета-(x,y:Вещ):Вещ;
(* Бета(x,y) = Г(x)*Г(y)/Г(x+y) *)
ПЕР
Lx,Ly,Lxy:Вещ;
SgnBeta:ЦЕЛ;
УКАЗ
Матем.код:=0;
SgnBeta:=ЗнакГ(x)*ЗнакГ(y)*ЗнакГ(x+y);
Lxy:=lnГ(x+y);
ЕСЛИ Матем.код # 0 ТО
ВОЗВРАТ 0
КОН;
Lx:=lnГ(x);
ЕСЛИ Матем.код # 0 ТО
ВОЗВРАТ SgnBeta*Матем.МАКСВЕЩ
КОН;
Ly:=lnГ(y);
ЕСЛИ Матем.код # 0 ТО
ВОЗВРАТ SgnBeta*Матем.МАКСВЕЩ
КОН;
ВОЗВРАТ SgnBeta*Матем.exp(Lx+Ly-Lxy)
КОН Бета;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА PSeries(a,b,x:Вещ):Вещ;
(* Степенной ряд для неполного бета интеграла.
* Используется когда b*x мало. *)
ПЕР
s,T,U,V,T1,z,Ai:Вещ;
n:ЦЕЛ;
УКАЗ
Ai:=1/a;
U:=(1-b)*x;
V:=U/(a+1);
T1:=V;
T:=U;
n:=2;
s:=0;
z:=Матем.МАШЕПС*Ai;
ПОКА МОДУЛЬ(V) > z ВЫП
U:=(n-b)*x/n;
T:=T*U;
V:=T/(a+n);
s:=s+V;
n:=n+1;
КОН;
s:=s+T1;
s:=s+Ai;
U:=a*Матем.ln(x);
ЕСЛИ (a+b < МАКСГ) И (МОДУЛЬ(U) < Матем.МАКСЛН) ТО
T:=Г(a+b)/(Г(a)*Г(b));
s:=s*T*Матем.степВВ(x,a)
ИНАЧЕ
T:=lnГ(a+b)-lnГ(a)-lnГ(b)+U+Матем.ln(s);
ЕСЛИ T < Матем.МИНЛН ТО
s:=0
ИНАЧЕ
s:=Матем.exp(T)
КОН
КОН;
ВОЗВРАТ s
КОН PSeries;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА CFrac1(a,b,x:Вещ):Вещ;
(* Цепная дробь для расширения # 1 в неполном бета интеграле. *)
ПЕР
Xk,Pk,Pkm1,Pkm2,Qk,Qkm1,Qkm2,K1,K2,K3,K4,K5,K6,K7,K8,r,T,ответ,Thresh:Вещ;
n:ЦЕЛ;
УКАЗ
K1:=a;
K2:=a+b;
K3:=a;
K4:=a+1;
K5:=1;
K6:=b-1;
K7:=K4;
K8:=a+2;
Pkm2:=0;
Qkm2:=1;
Pkm1:=1;
Qkm1:=1;
ответ:=1;
r:=1;
n:=0;
Thresh:=3*Матем.МАШЕПС;
ПОВТОРЯТЬ
Xk:=-(x*K1*K2)/(K3*K4);
Pk:=Pkm1+Pkm2*Xk;
Qk:=Qkm1+Qkm2*Xk;
Pkm2:=Pkm1;
Pkm1:=Pk;
Qkm2:=Qkm1;
Qkm1:=Qk;
Xk:=(x*K5*K6)/(K7*K8);
Pk:=Pkm1+Pkm2*Xk;
Qk:=Qkm1+Qkm2*Xk;
Pkm2:=Pkm1;
Pkm1:=Pk;
Qkm2:=Qkm1;
Qkm1:=Qk;
ЕСЛИ Qk # 0 ТО r:=Pk/Qk КОН;
ЕСЛИ r # 0 ТО
T:=МОДУЛЬ((ответ-r)/r);
ответ:=r
ИНАЧЕ
T:=1
КОН;
ЕСЛИ T < Thresh ТО
ВОЗВРАТ ответ;
КОН;
K1:=K1+1;
K2:=K2+1;
K3:=K3+2;
K4:=K4+2;
K5:=K5+1;
K6:=K6-1;
K7:=K7+2;
K8:=K8+2;
ЕСЛИ МОДУЛЬ(Qk)+МОДУЛЬ(Pk) > BIG ТО
Pkm2:=Pkm2*BIGINV;
Pkm1:=Pkm1*BIGINV;
Qkm2:=Qkm2*BIGINV;
Qkm1:=Qkm1*BIGINV
КОН;
ЕСЛИ (МОДУЛЬ(Qk) < BIGINV) ИЛИ (МОДУЛЬ(Pk) < BIGINV) ТО
Pkm2:=Pkm2*BIG;
Pkm1:=Pkm1*BIG;
Qkm2:=Qkm2*BIG;
Qkm1:=Qkm1*BIG
КОН;
n:=n+1
ДО n > 400;
Матем.код:=Матем.ОБНУЛЕНИЕ;
ВОЗВРАТ ответ
КОН CFrac1;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА CFrac2(a,b,x:Вещ):Вещ;
(* Цепная дробь для расширения # 2 в неполном бета интеграле. *)
ПЕР
Xk,Pk,Pkm1,Pkm2,Qk,Qkm1,Qkm2,K1,K2,K3,K4,K5,K6,K7,K8,r,T,z,ответ,Thresh:Вещ;
n:ЦЕЛ;
УКАЗ
K1:=a;
K2:=b-1;
K3:=a;
K4:=a+1;
K5:=1;
K6:=a+b;
K7:=a+1;
K8:=a+2;
Pkm2:=0;
Qkm2:=1;
Pkm1:=1;
Qkm1:=1;
z:=x/(1-x);
ответ:=1;
r:=1;
n:=0;
Thresh:=3*Матем.МАШЕПС;
ПОВТОРЯТЬ
Xk:=-(z*K1*K2)/(K3*K4);
Pk:=Pkm1+Pkm2*Xk;
Qk:=Qkm1+Qkm2*Xk;
Pkm2:=Pkm1;
Pkm1:=Pk;
Qkm2:=Qkm1;
Qkm1:=Qk;
Xk:=(z*K5*K6)/(K7*K8);
Pk:=Pkm1+Pkm2*Xk;
Qk:=Qkm1+Qkm2*Xk;
Pkm2:=Pkm1;
Pkm1:=Pk;
Qkm2:=Qkm1;
Qkm1:=Qk;
ЕСЛИ Qk # 0 ТО r:=Pk/Qk КОН;
ЕСЛИ r # 0 ТО
T:=МОДУЛЬ((ответ-r)/r);
ответ:=r
ИНАЧЕ
T:=1
КОН;
ЕСЛИ T < Thresh ТО
ВОЗВРАТ ответ
КОН;
K1:=K1+1;
K2:=K2-1;
K3:=K3+2;
K4:=K4+2;
K5:=K5+1;
K6:=K6+1;
K7:=K7+2;
K8:=K8+2;
ЕСЛИ МОДУЛЬ(Qk)+МОДУЛЬ(Pk) > BIG ТО
Pkm2:=Pkm2*BIGINV;
Pkm1:=Pkm1*BIGINV;
Qkm2:=Qkm2*BIGINV;
Qkm1:=Qkm1*BIGINV
КОН;
ЕСЛИ (МОДУЛЬ(Qk) < BIGINV) ИЛИ (МОДУЛЬ(Pk) < BIGINV) ТО
Pkm2:=Pkm2*BIG;
Pkm1:=Pkm1*BIG;
Qkm2:=Qkm2*BIG;
Qkm1:=Qkm1*BIG
КОН;
n:=n+1
ДО n > 400;
Матем.код:=Матем.ОБНУЛЕНИЕ;
ВОЗВРАТ ответ
КОН CFrac2;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА НепБета-(a,b,x:Вещ):Вещ;
ПЕР
A1,B1,x1,T,W,Xc,y:Вещ;
Flag:КЛЮЧ;
УКАЗ
Матем.код:=0;
ЕСЛИ (a <= 0) ИЛИ (b <= 0) ИЛИ (x < 0) ИЛИ (x > 1) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ (x = 0) ИЛИ (x = 1) ТО
ВОЗВРАТ x
КОН;
Flag:=ОТКЛ;
ЕСЛИ (b*x <= 1) И (x <= 0.95) ТО
T:=PSeries(a,b,x);
ЕСЛИ Flag ТО
ЕСЛИ T <= Матем.МАШЕПС ТО
T:=1-Матем.МАШЕПС
ИНАЧЕ
T:=1-T
КОН
КОН;
ВОЗВРАТ T
КОН;
W:=1-x;
(* Поменять местами < a > с < b >, если < x > больше среднего. *)
ЕСЛИ x > a/(a+b) ТО
Flag:=ВКЛ;
A1:=b;
B1:=a;
Xc:=x;
x1:=W
ИНАЧЕ
A1:=a;
B1:=b;
Xc:=W;
x1:=x
КОН;
ЕСЛИ Flag И (B1*x1 <= 1) И (x1 <= 0.95) ТО
T:=PSeries(A1,B1,x1);
ЕСЛИ Flag ТО
ЕСЛИ T <= Матем.МАШЕПС ТО
T:=1-Матем.МАШЕПС
ИНАЧЕ
T:=1-T
КОН
КОН;
ВОЗВРАТ T
КОН;
(* Применить расширение для лучшей сходимости *)
y:=x1*(A1+B1-2)-(A1-1);
ЕСЛИ y < 0 ТО
W:=CFrac1(A1,B1,x1)
ИНАЧЕ
W:=CFrac2(A1,B1,x1)/Xc
КОН;
(* Умножение w на
* _ _
* xa(1-x)bГ(a+b)/( a | (a) | (b) ) *)
y:=A1*Матем.ln(x1);
T:=B1*Матем.ln(Xc);
ЕСЛИ (A1+B1 < МАКСГ) И (МОДУЛЬ(y) < Матем.МАКСЛН) И (МОДУЛЬ(T) < Матем.МАКСЛН) ТО
T:=Матем.степВВ(Xc,B1);
T:=T*Матем.степВВ(x1,A1);
T:=T/A1;
T:=T*W;
T:=T*Г(A1+B1)/(Г(A1)*Г(B1));
ИНАЧЕ
(* Resort to logarithms *)
y:=y+T+lnГ(A1+B1)-lnГ(A1)-lnГ(B1)+Матем.ln(W/A1);
ЕСЛИ y < Матем.МИНЛН ТО
T:=0
ИНАЧЕ
T:=Матем.exp(y)
КОН
КОН;
ЕСЛИ Flag ТО
ЕСЛИ T <= Матем.МАШЕПС ТО
T:=1-Матем.МАШЕПС
ИНАЧЕ
T:=1-T
КОН
КОН;
ВОЗВРАТ T;
КОН НепБета;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА erf-(x:Вещ):Вещ;
УКАЗ
ЕСЛИ x < 0 ТО
ВОЗВРАТ -НепГ(0.5,Матем.кв(x))
КОН;
ВОЗВРАТ НепГ(0.5,Матем.кв(x))
КОН erf;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА erfc-(x:Вещ):Вещ;
УКАЗ
ЕСЛИ x < 0 ТО
ВОЗВРАТ 1+НепГ(0.5,Матем.кв(x))
КОН;
ВОЗВРАТ ДопГ(0.5,Матем.кв(x))
КОН erfc;
(*----------------------------------------------------------------------------*
* Распределение вероятностей
*----------------------------------------------------------------------------*)
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА БиномВ-(n:ЦЕЛ; p:Вещ; k:ЦЕЛ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (p < 0) ИЛИ (p > 1) ИЛИ (n <= 0) ИЛИ (n < k) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ k = 0 ТО
ВОЗВРАТ Матем.степВЦ(1-p,n)
КОН;
ЕСЛИ k = n ТО
ВОЗВРАТ Матем.степВЦ(p,n)
КОН;
ВОЗВРАТ БиномКоэф(n,k)*Матем.степВЦ(p,k)*Матем.степВЦ(1-p,n-k)
КОН БиномВ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА БиномП-(n:ЦЕЛ; p:Вещ; k:ЦЕЛ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (p < 0) ИЛИ (p > 1) ИЛИ (n <= 0) ИЛИ (n < k) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ k = 0 ТО
ВОЗВРАТ Матем.степВЦ(1-p,n)
КОН;
ЕСЛИ k = n ТО
ВОЗВРАТ 1
КОН;
ВОЗВРАТ 1-НепБета(k+1,n-k,p)
КОН БиномП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ПуассонВ-(Mu:Вещ; k:ЦЕЛ):Вещ;
ПЕР
p:Вещ;
i:ЦЕЛ;
УКАЗ
Матем.код:=0;
ЕСЛИ (Mu <= 0) ИЛИ (k < 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ k = 0 ТО
ВОЗВРАТ Матем.exp(-Mu)
КОН;
p:=Mu;
ОТ i:=2 ДО k ВЫП
p:=p*Mu/i
КОН;
ВОЗВРАТ Матем.exp(-Mu)*p
КОН ПуассонВ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ПуассонП-(Mu:Вещ; k:ЦЕЛ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (Mu <= 0) ИЛИ (k < 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ k = 0 ТО
ВОЗВРАТ Матем.exp(-Mu)
КОН;
ВОЗВРАТ ДопГ(k+1,Mu)
КОН ПуассонП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА НормПл-(x:Вещ):Вещ;
УКАЗ
ВОЗВРАТ Матем.exp(-0.5*Матем.кв(x))/Матем.КВКОР2ПИ
КОН НормПл;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА НормП-(x:Вещ):Вещ;
УКАЗ
ВОЗВРАТ 0.5*(1+erf(x*Матем.КВКОР2/2))
КОН НормП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ОбрНорм-(p:Вещ):Вещ;
(* Обратная функция нормального распределения
* Возвращает такое < x >, при котором площадь под графиком функции плотности
* вероятности Гауссова распределения (от -БСК до < x >) равна < p >. *)
ПЕР
x,y,z,y2,x0,x1:Вещ;
Code:ЦЕЛ;
УКАЗ
ЕСЛИ (p <= 0) ИЛИ (p >= 1) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0;
КОН;
Code:=1;
y:=p;
ЕСЛИ y > (1-0.13533528323661269189) ТО (* 0.135... = Матем.exp(-2) *)
y:=1-y;
Code:=0;
КОН;
ЕСЛИ y > 0.13533528323661269189 ТО
y:=y-0.5;
y2:=y*y;
x:=y+y*(y2*Горнер0(y2,INP0,7)/Горнер1(y2,INQ0,7));
x:=x*Матем.КВКОР2ПИ;
ВОЗВРАТ x
КОН;
x:=Матем.квкор(-2*Матем.ln(y));
x0:=x-Матем.ln(x)/x;
z:=1/x;
ЕСЛИ x < 8 ТО
x1:=z*Горнер0(z,INP1,9)/Горнер1(z,INQ1,9)
АЕСЛИ x < 32 ТО
x1:=z*Горнер0(z,INP2,7)/Горнер1(z,INQ2,7)
ИНАЧЕ
x1:=z*Горнер0(z,INP3,7)/Горнер1(z,INQ3,7)
КОН;
x:=x0-x1;
ЕСЛИ Code # 0 ТО x:=-x КОН;
ВОЗВРАТ x
КОН ОбрНорм;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА НормВ-(x:Вещ):Вещ;
ПЕР
a:Вещ;
УКАЗ
a:=МОДУЛЬ(x);
Матем.код:=0;
ЕСЛИ a = 0 ТО
ВОЗВРАТ 1
КОН;
ЕСЛИ a < 1 ТО
ВОЗВРАТ 1-erf(a*Матем.КВКОР2/2)
КОН;
ВОЗВРАТ erfc(a*Матем.КВКОР2/2)
КОН НормВ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА СтьюдентПл-(Nu:ЦЕЛ; x:Вещ):Вещ;
ПЕР
l,p,q:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ Nu < 1 ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
p:=0.5*(Nu+1);
q:=0.5*Nu;
l:=lnГ(p)-lnГ(q)-0.5*Матем.ln(Nu*Матем.ПИ)-p*Матем.ln(1+Матем.кв(x)/Nu);
ВОЗВРАТ Матем.exp(l)
КОН СтьюдентПл;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА СтьюдентП-(Nu:ЦЕЛ; x:Вещ):Вещ;
ПЕР
f:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ Nu < 1 ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ x = 0 ТО
ВОЗВРАТ 0.5
КОН;
f:=0.5*НепБета(0.5*Nu,0.5,Nu/(Nu+Матем.кв(x)));
ЕСЛИ x < 0 ТО
ВОЗВРАТ f
ИНАЧЕ
ВОЗВРАТ 1-f
КОН
КОН СтьюдентП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА СтьюдентВ-(Nu:ЦЕЛ; x:Вещ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ Nu < 1 ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ НепБета(0.5*Nu,0.5,Nu/(Nu+Матем.кв(x)))
КОН СтьюдентВ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Хи2Пл-(Nu:ЦЕЛ; x:Вещ):Вещ;
УКАЗ
Матем.код:=0;
ВОЗВРАТ ГПл(0.5*Nu,0.5,x)
КОН Хи2Пл;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Хи2П-(Nu:ЦЕЛ; x:Вещ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (Nu < 1) ИЛИ (x <= 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ НепГ(0.5*Nu,0.5*x)
КОН Хи2П;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА Хи2В-(Nu:ЦЕЛ; x:Вещ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (Nu < 1) ИЛИ (x <= 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ ДопГ(0.5*Nu,0.5*x)
КОН Хи2В;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ФишерПл-(Nu1,Nu2:ЦЕЛ; x:Вещ):Вещ;
ПЕР
p1,p2,r,s,l:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (Nu1 < 1) ИЛИ (Nu2 < 1) ИЛИ (x <= 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
r:=Nu1 ДЕЛИТЬ Nu2;
p1:=0.5*Nu1;
p2:=0.5*Nu2;
s:=p1+p2;
l:=lnГ(s)-lnГ(p1)-lnГ(p2)+p1*Матем.ln(r);
l:=l+(p1-1)*Матем.ln(x)-s*Матем.ln(1+r*x);
ВОЗВРАТ Матем.exp(l)
КОН ФишерПл;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ФишерП-(Nu1,Nu2:ЦЕЛ; x:Вещ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (Nu1 < 1) ИЛИ (Nu2 < 1) ИЛИ (x <= 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ 1-НепБета(0.5*Nu2,0.5*Nu1,Nu2/(Nu2+Nu1*x))
КОН ФишерП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ФишерВ-(Nu1,Nu2:ЦЕЛ; x:Вещ):Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (Nu1 < 1) ИЛИ (Nu2 < 1) ИЛИ (x <= 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ НепБета(0.5*Nu2,0.5*Nu1,Nu2/(Nu2+Nu1*x))
КОН ФишерВ;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ЭкспПл-(a,x:Вещ):Вещ;
УКАЗ
ЕСЛИ (a <= 0) ИЛИ (x < 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ a*Матем.exp(-a*x)
КОН ЭкспПл;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ЭкспП-(a,x:Вещ):Вещ;
УКАЗ
ЕСЛИ (a <= 0) ИЛИ (x < 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ 1-Матем.exp(-a*x)
КОН ЭкспП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА БетаПл-(a,b,x:Вещ):Вещ;
ПЕР
l:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (a <= 0) ИЛИ (b <= 0) ИЛИ (x < 0) ИЛИ (x > 1) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ x = 0 ТО
ЕСЛИ a < 1 ТО
Матем.код:=Матем.ОСОБЕННОСТЬ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ 0
КОН;
ЕСЛИ x = 1 ТО
ЕСЛИ b < 1 ТО
Матем.код:=Матем.ОСОБЕННОСТЬ;
ВОЗВРАТ 0
КОН;
ВОЗВРАТ 0
КОН;
l:=lnГ(a+b)-lnГ(a)-lnГ(b);
l:=l+(a-1)*Матем.ln(x)+(b-1)*Матем.ln(1-x);
ВОЗВРАТ Матем.exp(l)
КОН БетаПл;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА БетаП-(a,b,x:Вещ):Вещ;
УКАЗ
ВОЗВРАТ НепБета(a,b,x)
КОН БетаП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ГПл-(a,b,x:Вещ):Вещ;
ПЕР
l:Вещ;
УКАЗ
Матем.код:=0;
ЕСЛИ (a <= 0) ИЛИ (b <= 0) ИЛИ (x < 0) ТО
Матем.код:=Матем.ВНЕОБЛАСТИ;
ВОЗВРАТ 0
КОН;
ЕСЛИ x = 0 ТО
ЕСЛИ a < 1 ТО
Матем.код:=Матем.ОСОБЕННОСТЬ;
ВОЗВРАТ 0
КОН;
ЕСЛИ a = 1 ТО
ВОЗВРАТ b
КОН;
ВОЗВРАТ 0
КОН;
l:=a*Матем.ln(b)-lnГ(a)+(a-1)*Матем.ln(x)-b*x;
ВОЗВРАТ Матем.exp(l)
КОН ГПл;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
ЗАДАЧА ГП-(a,b,x:Вещ):Вещ;
УКАЗ
ВОЗВРАТ НепГ(a,b*x)
КОН ГП;
(*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*)
УКАЗ
ТаблФакт[0]:=1;
ТаблФакт[1]:=1;
ТаблФакт[2]:=2;
ОТ i:=3 ДО ТАБЛФАКТ-1 ВЫП
ТаблФакт[i]:=ТаблФакт[i-1]*i
КОН;
STIR[0]:= 7.147391378143610789273D-4;
STIR[1]:=-2.363848809501759061727D-5;
STIR[2]:=-5.950237554056330156018D-4;
STIR[3]:= 6.989332260623193171870D-5;
STIR[4]:= 7.840334842744753003862D-4;
STIR[5]:=-2.294719747873185405699D-4;
STIR[6]:=-2.681327161876304418288D-3;
STIR[7]:= 3.472222222230075327854D-3;
STIR[8]:= 8.333333333333331800504D-2;
STIR[9]:= 0;
GAMS[0]:=-1.193945051381510095614D-3;
GAMS[1]:= 7.220599478036909672331D-3;
GAMS[2]:=-9.622023360406271645744D-3;
GAMS[3]:=-4.219773360705915470089D-2;
GAMS[4]:= 1.665386113720805206758D-1;
GAMS[5]:=-4.200263503403344054473D-2;
GAMS[6]:=-6.558780715202540684668D-1;
GAMS[7]:= 5.772156649015328608253D-1;
GAMS[8]:= 1.000000000000000000000D0;
GAMS[9]:= 0;
GAMSN[0]:= 1.133374167243894382010D-3;
GAMSN[1]:= 7.220837261893170325704D-3;
GAMSN[2]:= 9.621911155035976733706D-3;
GAMSN[3]:=-4.219773343731191721664D-2;
GAMSN[4]:=-1.665386113944413519335D-1;
GAMSN[5]:=-4.200263503402112910504D-2;
GAMSN[6]:= 6.558780715202536547116D-1;
GAMSN[7]:= 5.772156649015328608727D-1;
GAMSN[8]:=-1.000000000000000000000D0;
GAMSN[9]:= 0;
STIRP[0]:= 4.885026142432270781165D-3;
STIRP[1]:=-1.880801938119376907179D-3;
STIRP[2]:= 8.412723297322498080632D-4;
STIRP[3]:=-5.952345851765688514613D-4;
STIRP[4]:= 7.936507795855070755671D-4;
STIRP[5]:=-2.777777777750349603440D-3;
STIRP[6]:= 8.333333333333331447505D-2;
STIRP[7]:= 0;
GAMP[0]:= 4.212760487471622013093D-5;
GAMP[1]:= 4.542931960608009155600D-4;
GAMP[2]:= 4.092666828394035500949D-3;
GAMP[3]:= 2.385363243461108252554D-2;
GAMP[4]:= 1.113062816019361559013D-1;
GAMP[5]:= 3.629515436640239168939D-1;
GAMP[6]:= 8.378004301573126728826D-1;
GAMP[7]:= 1.000000000000000000009D0;
GAMP[8]:= 0;
GAMQ[0]:=-1.397148517476170440917D-5;
GAMQ[1]:= 2.346584059160635244282D-4;
GAMQ[2]:=-1.237799246653152231188D-3;
GAMQ[3]:=-7.955933682494738320586D-4;
GAMQ[4]:= 2.773706565840072979165D-2;
GAMQ[5]:=-4.633887671244534213831D-2;
GAMQ[6]:=-2.243510905670329164562D-1;
GAMQ[7]:= 4.150160950588455434583D-1;
GAMQ[8]:= 9.999999999999999999908D-1;
GAMQ[9]:= 0;
LNGAMP[0]:=-2.163690827643812857640D3;
LNGAMP[1]:=-8.723871522843511459790D4;
LNGAMP[2]:=-1.104326814691464261197D6;
LNGAMP[3]:=-6.111225012005214299996D6;
LNGAMP[4]:=-1.625568062543700591014D7;
LNGAMP[5]:=-2.003937418103815175475D7;
LNGAMP[6]:=-8.875666783650703802159D6;
LNGAMP[7]:= 0;
LNGAMP[8]:= 0;
LNGAMP[9]:= 0;
LNGAMQ[0]:=-5.139481484435370143617D2;
LNGAMQ[1]:=-3.403570840534304670537D4;
LNGAMQ[2]:=-6.227441164066219501697D5;
LNGAMQ[3]:=-4.814940379411882186630D6;
LNGAMQ[4]:=-1.785433287045078156959D7;
LNGAMQ[5]:=-3.138646407656182662088D7;
LNGAMQ[6]:=-2.099336717757895876142D7;
LNGAMQ[7]:= 0;
LNGAMQ[8]:= 0;
LNGAMQ[9]:= 0;
INP0[0]:= 8.779679420055069160496D-3;
INP0[1]:=-7.649544967784380691785D-1;
INP0[2]:= 2.971493676711545292135D0;
INP0[3]:=-4.144980036933753828858D0;
INP0[4]:= 2.765359913000830285937D0;
INP0[5]:=-9.570456817794268907847D-1;
INP0[6]:= 1.659219375097958322098D-1;
INP0[7]:=-1.140013969885358273307D-2;
INP0[8]:= 0;
INP0[9]:= 0;
INQ0[0]:=-5.303846964603721860329D0;
INQ0[1]:= 9.908875375256718220854D0;
INQ0[2]:=-9.031318655459381388888D0;
INQ0[3]:= 4.496118508523213950686D0;
INQ0[4]:=-1.250016921424819972516D0;
INQ0[5]:= 1.823840725000038842075D-1;
INQ0[6]:=-1.088633151006419263153D-2;
INQ0[7]:= 0;
INQ0[8]:= 0;
INQ0[9]:= 0;
INP1[0]:= 4.302849750435552180717D0;
INP1[1]:= 4.360209451837096682600D1;
INP1[2]:= 9.454613328844768318162D1;
INP1[3]:= 9.336735653151873871756D1;
INP1[4]:= 5.305046472191852391737D1;
INP1[5]:= 1.775851836288460008093D1;
INP1[6]:= 3.640308340137013109859D0;
INP1[7]:= 3.691354900171224122390D-1;
INP1[8]:= 1.403530274998072987187D-2;
INP1[9]:= 1.377145111380960566197D-4;
INQ1[0]:= 2.001425109170530136741D1;
INQ1[1]:= 7.079893963891488254284D1;
INQ1[2]:= 8.033277265194672063478D1;
INQ1[3]:= 5.034715121553662712917D1;
INQ1[4]:= 1.779820137342627204153D1;
INQ1[5]:= 3.845554944954699547539D0;
INQ1[6]:= 3.993627390181238962857D-1;
INQ1[7]:= 1.526870689522191191380D-2;
INQ1[8]:= 1.498700676286675466900D-4;
INQ1[9]:= 0;
INP2[0]:= 3.244525725312906932464D0;
INP2[1]:= 6.856256488128415760904D0;
INP2[2]:= 3.765479340423144482796D0;
INP2[3]:= 1.240893301734538935324D0;
INP2[4]:= 1.740282292791367834724D-1;
INP2[5]:= 9.082834200993107441750D-3;
INP2[6]:= 1.617870121822776093899D-4;
INP2[7]:= 7.377405643054504178605D-7;
INP2[8]:= 0;
INP2[9]:= 0;
INQ2[0]:= 6.021509481727510630722D0;
INQ2[1]:= 3.528463857156936773982D0;
INQ2[2]:= 1.289185315656302878699D0;
INQ2[3]:= 1.874290142615703609510D-1;
INQ2[4]:= 9.867655920899636109122D-3;
INQ2[5]:= 1.760452434084258930442D-4;
INQ2[6]:= 8.028288500688538331773D-7;
INQ2[7]:= 0;
INQ2[8]:= 0;
INQ2[9]:= 0;
INP3[0]:= 2.020331091302772535752D0;
INP3[1]:= 2.133020661587413053144D0;
INP3[2]:= 2.114822217898707063183D-1;
INP3[3]:=-6.500909615246067985872D-3;
INP3[4]:=-7.279315200737344309241D-4;
INP3[5]:=-1.275404675610280787619D-5;
INP3[6]:=-6.433966387613344714022D-8;
INP3[7]:=-7.772828380948163386917D-11;
INP3[8]:= 0;
INP3[9]:= 0;
INQ3[0]:= 2.278210997153449199574D0;
INQ3[1]:= 2.345321838870438196534D-1;
INQ3[2]:=-6.916708899719964982855D-3;
INQ3[3]:=-7.908542088737858288849D-4;
INQ3[4]:=-1.387652389480217178984D-5;
INQ3[5]:=-7.001476867559193780666D-8;
INQ3[6]:=-8.458494263787680376729D-11;
INQ3[7]:= 0;
INQ3[8]:= 0;
INQ3[9]:= 0
КОН Функц.
|
|
Программируем по-русски